% Construct the finite difference matrix
N = 256;
ntotal = N^2;
epsilon = 1;


h = (1/N)/epsilon;
[x,y] = ndgrid( (0:(N-1))/N, (0:(N-1))/N );
% V = 5d-1*sin(2*pi / lc * (x/epsilon)) .* ...
  % sin(2*pi / lc * (y/epsilon) ) + 1d-3 * (rand(size(x))-0.5) .* ...
  % (rand(size(y))-0.5) ;
% V = zeros( N, N );
% V = 1d-1 * reshape(rand(ntotal,1), N, N);
V = 1 + rand(N,N);


% nx = N;
% ny = N;
% display('Generate a 2-d tight binding matrix');
% H = sparse(nx*nx, nx*nx);
% for i1 = 1 : nx
  % for j1 = 1 : ny
    % pos = nx*(i1-1) + j1;
    % pos1 = nx * (i1-1) + mod(j1,ny)+1;
    % pos2 = nx * (i1-1) + mod(j1-2,ny)+1;
    % pos3 = nx * (mod(i1,nx)) + j1;
    % pos4 = nx * (mod(i1-2,nx)) + j1;
    % H(pos, pos1) = -1;
    % H(pos, pos2) = -1;
    % H(pos, pos3) = -1;
    % H(pos, pos4) = -1;
    % H(pos, pos) = 4;
  % end
% end
% H = H / (h^2);
% for j1 = 1 : nx
 % for i1 = 1 : ny
   % pos = nx*(j1-1) + i1;
   % H(pos, pos) = H(pos, pos) + V(i1,j1);
 % end
% end

% display('Diagonalizing');
% NSite = (N/lc)^2;
% opts.disp = 0;
% [EF, D] = eigs( H, NSite+10, 'sm', opts );
% D = diag(D);
% SD = sort(D);

